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A two-electron one-dimensional model of a heteroatomic molecule composed of two open-shell 
atoms is considered. Including only two electrons isolates and examines the effect that the highest 
occupied molecular orbital has on the Kohn-Sham potential as the molecule dissociates. We repro- 
duce the characteristic step and peak that previous high-level wavefunction methods have shown to 
exist for real molecules in the low-density internuclear region. The simplicity of our model enables 
us to investigate in detail their development as a function of bond-length, with little computational 
effort, and derive properties of their features in the dissociation limit. We show that the onset of 
the step is coincident with the internuclear separation at which an avoided crossing between the 
ground-state and lowest charge-transfer excited state is approached. Although the step and peak 
features have little effect on the ground-state energetics, we discuss their important consequences 
for dynamics and response. 



I. INTRODUCTION 

The unprecedented balance between accuracy and 
efficiency of density functional theory (DFT)^"'' is 
due in large part to the discoveries of John Perdew. 
The mapping of the true system of interacting elec- 
trons to a fictitious one in which the electrons don't 
interact, yet reproduce the true electron density, re- 
quires accurate approximations for the exchange- 
correlation (xc) potential, which remained elusive 
until the developments in the 1980's of Perdew and 
co-workers. Understanding and incorporating exact 
conditions and physical principles underlie the ro- 
bustness and reliability of Perdew's functionals. In 
this spirit, we study here the structure of the exact 
xc potential as a molecule dissociates, whose land- 
scape of steps and peaks Perdew was one of the first 
to explore. 

In DPT, one solves self-consistently the Kohn- 
Sham (KS) single-particle equations 



according to 



1. 



(1) 



where Vs [p] (r) is the KS potential, a functional of the 
ground-state electronic density, p. (Atomic units, 
— Ti — me — I, are used throughout this pa- 
per). It is usually written as the sum Ws[p](r) = 
t^Gxt[p](r) + WH[p](r) + Uxc[p](r), where Uext(r) is the 
potential due to the nuclei, i^nCr) = / d^r' p{v')/\r — 
r'l is the classical Hartree potential and Wxc(r) is 
the exchange-correlation (xc) potential, incorporat- 
ing the remaining many-body effects in a one-body 
potential. The KS orbitals yield the true density 
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p(r)=5^|0,(r) 



(2) 



In principle, the ground-state density and all static 
properties of the true interacting system are exactly 
recovered, but in practice approximations are needed 
for the unknown xc potential Wa;c[p] as a functional of 
the density. Typically, semi- local functionals, such 
as GGA's^ and meta-GGA's'' give good energies and 
structural properties at equilibrium molecular ge- 
ometries; the non-empirical constructions of Perdew 
and co-workers impart a reliability to the descrip- 
tion of diverse systems and properties. However, 
GGA's do not perform so well for weakly-coupled 
sub-systems. Notably, recent work has been very 
successful in describing van der Waal's forces using 
sophisticated non-local approximations in DFT^"^-"^. 
For molecules dissociating into open-shells, the fail- 
ure of semi-local approximations becomes drastic, 
yielding unphysical fractional charges on the sep- 
arated species^^"^^. This problem was first high- 
lighted by Perdew^^'^'*, motivated by the observation 
of Slater^^ that his "Xa" method yields a similar re- 
sult. 

Figure 7 of Ref.^'* shows that the exact xc poten- 
tial develops a "region of positive constant" around 
the atom with the "tighter density distribution" , in 
the limit of infinite separation, using a simple one- 
dimensional model. More generally, the effect of 
molecular dissociation on the ground-state xc po- 
tential for the case of real diatomic closed shell 
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molecules consisting of open shell atoms has been 

studied systematically, by Bacrcnds, Gritscnko, and 
co-workers, in a scries of papcrs^®"^^. In Ref.^®, the 
simplest case of this, the two-electron H2 molecule 
was studied. The absence of long-range left-right 
correlation in Hartrec-Fock, renders its potential 
overly repulsive near the nuclei, leading to an overly 
diffuse density. A highly accurate xc potential was 
constructed from correlated CISD first- and second- 
order density matrices in Rcf.^®, and the resulting 
correlation potential was shown to considerably re- 
duce the repulsion at the nuclei. It was also shown 
that the xc potential develops a sharp maximum 
("peak") at the bond midpoint. A very thorough 
analysis of the KS potential in stretched H2 was 
performed later in Ref.^^, where the effect of dif- 
ferent approximate constructions for the KS orbital 
was investigated and explained in detail (see also 
Sec. IV B). Using an iterative method introduced in 
Ref.^^, Ref.^^ was the first to construct molecular 
KS potentials for more than two electrons from cor- 
related densities. They studied LiH (and H2) and 
found significant differences with the local density 
approximation (LDA) at large separations. Ref.-*^^ 
calculated the xc potential for the monohydrides XH 
(X=Li,B,F), analysing its structure via a decompo- 
sition, or "partitioning" of v^c into various "energy" 
and "response" components related to the electronic 
structure^®'^*'^°. It was shown that left-right corre- 
lation leads to a build-up in the xc potential around 
the H atom (a "step" , as was observed in the simple 
model of Ref.^'*). The "peak" present in the bond 
mid-point of H2 was found, in the case of the mono- 
hydrides, to shifted toward the H atom while becom- 
ing significantly smaller due to the presence of core 
electrons softening the left-right correlation effects. 
The partitioning scheme (reviewed in Sec. II), which 
had earlier bec;n used to examine atomic xc poten- 
tials^^'^", proved to be a particularly useful tool in 
the analysis, providing insight into the origin of the 
peak and step structures. 

Molecular dissociation in DFT is particularly rele- 
vant when considering time-dependent processes and 
nuclear dynamics on potential energy surfaces. The 
advent of time-dependent density functional theory 
(TDDFT) 23-25 allows for a density-functional de- 
scription of full electron dynamics and here accu- 
rate long-range potentials are an important ingre- 
dient for many applications, eg. photo-dissociation 
dynamics, excitations in large molecules, including 
charge-transfer, and molecular transport. A recent 
paper^^ discussed promising aspects as well as chal- 
lenges in getting accurate excited energy surfaces 
from TDDFT; certainly it is important to get the 
ground-state potential energy surface correct. 

In the present paper, we study the xc potential of a 



dissociating closed-shell hetero-atomic molecule con- 
sisting of two open-shell atoms by analysing a simple 
one-dimensional model of two different "one-electron 
atoms" . The two "electrons" and two "nuclei" inter- 
act via soft-Coulomb interactions with the softening 
parameters chosen to approximate certain proper- 
ties of the real LiH molecule. This simple model 
allows numerically exact solution at a wide range 
of separations with little computational effort, while 
reproducing the essential features, from the point of 
view of molecular dissociation, of the xc potential for 
real three-dimensional molecules. It allows some an- 
alytic treatment of these features that yields further 
insight into the "step" and "peak" structures men- 
tioned above; for example, predicting the asymptotic 
height and position of the peak and an explanation 
of why such a structure, that hardly affects the en- 
ergetics, must be there. A detailed examination of 
the stretched bond-length where the step begins to 
appear, reveals a correlation with the position of the 
avoided crossing between the ground-state and low- 
est charge-transfer states. We explain why. 

A two-electron model isolates the effects due to 
the valence electrons, which play the major role in 
dissociative processes, without additional potential- 
features arising from core electrons. In the KS de- 
scription of dissociation, the major role is played by 
the KS HOMO, which, in the case of open-shell frag- 
ments, is delocalized across the molecule. By includ- 
ing only the HOMO in our model, we isolate and 
examine effects on the dissociating potential energy 
surface due solely to this most important orbital. 
The model is presented in Sec. HI while Sec. IV con- 
tain the numerical and analytic results. 

We may draw conclusions from this simple model 

about real three-dimensional molecules composed 
of open-shell fragments of general odd electron- 
number, but with a little caution: we shall find 
quantitative differences due to the lack of core elec- 
trons in our model, and the much longer effective 
range of the soft-Coulomb interaction in ID com- 
pared to the true 3D Coulomb interaction. The 
soft-Coulomb interaction is used in many inter- 
esting investigations of strong- field dynamics2'''-32 ^ 
and recently, in the context of TDDFT^^.s^. ^i^gge 
models capture the essence of phenomena such as 
non-sequential double-ionization, and laser-induced 
electron-recollision. The peak and step are challeng- 
ing features for approximations to capture, and are 
lacking in almost all functionals used today. Being 
in a low-density region, the peak structure has neg- 
ligible energetic consequence, however it does play 
a role when response or full dynamics is considered: 
for example, it reduces the (hyper-)polarizability of 
long-chain systems"^"''. As they represent barriers to 
electron transport, the work here is also relevant 
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to one-dimensional transport calculations in molec- 
ular wires^^, although this fact has not been dis- 
cussed before, so perhaps not yet been fully appre- 
ciated. The step structure, essential to avoid the 
fractional charge problem, has severe consequences 
for the structure of the TDDFT xc kernel as we shall 
discuss in Section V. 



II. DECOMPOSITION OF THE XC 
POTENTIAL 

The partitioning of the xc potential^^"^" was mo- 
tivated by first decomposing the xc energy compo- 
nents into "potential" and "kinetic" terms of the 
form: 

SxcW = Wxo[p\+To[p\ where 
W.o[p] = \j p(r)t;^f [p](r)dV and 

Tc[p] = J pirV^PKr) (3) 

implicitly define "hole" and "correlation kinetic" po- 
tentials, ^xc'' and w^'". The total xc potential is then 
partitioned into three components: 

^;xc [p] (r) = v^t [P] (r) + t^ckin [p] (r) + t;.esp [p] (r) (4) 



..hole 



[p](r) is the Coulomb potential of the xc hole: 

^hole(r) = J '-^^^dh, (5) 



where the xc hole Pxc(ri,r2) is defined through 
the pair density P(ri,r2) (joint probability of find- 
ing an electron at ri while another is at r2), via 
P(ri,r2) = p{ri){p{r2) + Pxc(ri, r2)). When added 
to the Hartree potential, t;xc''(r) + *^H(r) represents 
the average repulsion an electron at position r expe- 
riences due to the other (N-1) electrons in the sys- 
tem. In terms of the conditional probability ampli- 
tude, whose square gives the probability of finding 
the other (N-1) electrons in the system with space- 
spin coordinates X2, X3, x^v when an electron is 
known to be at position ri : 

I ^ *(xi,X2,...,Xiv) , , 
$(si,X2,...,XAr|ri) = ^ j= (6) 



V N 



where \I'(xi,X2, ...,XAr) is the interacting many- 
electron wavefunction, we have 

v^t{r)+vn{r)^ J $*(si,x2,...,x^|r) 



^ 1 

5^ It- _ 



.1=2 



$(s, X2, Xjv|r)dsidx2...rfx^) 



The second term in Eq. (4), Wc,km[/9](r), is the cor- 
relation contribution to the kinetic component of the 
xc potential. It is the difference of the kinetic com- 
ponents of the interacting and non-interacting KS 
systems: 



Wc,kiii(r) = Wkiii(r) - t;s,kin(r) 



(8) 



where the kinetic components may be written in 
terms of the conditional probability amplitude: 

t^kin(ri) = |Vi$(si,X2,...,XAr|ri)|^rfsidx2...dxAr 
Vi-Vip(r'i,ri)|r,=ri [Vp(r0]2 



P(ri) 



8p(ri)2 



(9) 



and 

^^s,kin(ri) 



^ / |Vi$s(si,X2,...,Xjv|ri)|^dsidx2...dxAr 



1 ^ 



'/'^(I■l) 



(10) 



In Eq. (9) p{r[,ri) is the first-order spin- 
summed reduced density-matrix, and in Eq. (10) 
$s(si,X2, ...,xjv|ri) is the conditional probability 
amplitude of the KS system, which is defined as in 
Eq. 6 but with the KS single Slater determinant, 
whose orbitals are ^i{r), replacing the full many- 
electron wavefunction. 

The final term in Eq. (4) is the so-called response 
potential. It may be further partitioned into terms 
representing the response of the xc hole, and the re- 
sponse of the correlation kinetic potential i^'^^, but 
we will not pursue this further decomposition here. 

For two electrons in a spin-singlet, many of these 
expressions simplify considerably, as the KS single- 
Slater determinant consists of one doubly-occupied 
spatial orbital, (/'o(r)) which is equal to square root 
of half the density: 



(11) 



Substituting into Eq. 1, we can solve explicitly for 
the KS potential as a functional of the density: 



(12) 



where / is the first ionization potential of the system. 
From Eq. (10), i's,kin(r) = and Eq. 8 reduces to: 



^^c,kin(r) = t;kin(r) 



(13) 



Also, for two electrons, ?;x(r) = Wx'^°'''(r) = 
— ^Vnir). The exchange component to the response 
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potential, Uresp is zero. We may therefore write: 

v^{r) = -v^{r)/2 = v'^°'%r) (14) 
^;c(r) = vJ;°'^(r)+t;c,ki„(r)+Vc,resp(r) (15) 

As was found in Refs.^*'^^, as a heteroatomic 
molecule dissociates, a step structure in the low- 
density bond midpoint region arises in the response 
component Wresp, accompanied by a peak structure 
in the kinetic component Uc,kin of the xc potential. 
(See also Figures 3 and 4). 

Consider now a simplified description of the 
molecule that has includes just one electron on each 
atom. As explained in Ref.^^, and reflected in Eq. 9, 
^'kin(r) depends on the gradient of the conditional 
probability amplitude, so describes how strongly 
the motion of an electron at reference position r 
is correlated with the other electrons in the sys- 
tem. For r near one of the nuclei, the reference 
electron moves in a potential dominated by the nu- 
clear potential, and the conditional amplitude $ re- 
duces to the atomic HOMO of the other atom, and 
doesn't change for small changes around r; hence 
Vhin goes to zero. But in the internuclear region, 
the motion of the two electrons becomes correlated: 
as the reference position moves from one nucleus to 
the other, the conditional probability of finding the 
other switches from being towards one atom to the 
other, and so Ukin peaks. 

The origin of the step structure was also analyzed 
extensively in^^ and shown to arise in the correla- 
tion component of the response potential, v^c^^. It 
was discovered earlier^^'^^, in relation to the deriva- 
tive discontinuity^"^'^^'"'^''*", that the correlation po- 
tential for a long-range molecule composed of two 
open-shell atoms must have a step in the midpoint 
region, whose size is such that the atomic HOMO 
orbital energies re-align. From Koopman's theorem, 
the HOMO energy equals the ionization potential; 
therefore the step has a size AI = I2 — h where 
72,1 is the larger(smaller) ionization potential of the 
two atoms, raising the potential of the more tightly 
bound atom. Far away from the molecule the po- 
tential near this atom steps back down to zero. A 
simple way to understand the origin of the step is 
to realise that had the step not been there, then 
one could lower the ground-state energy of the long- 
range molecule by transferring a fraction of charge 
from the atom with the higher ionization potential 
to that with the lower, leading to the molecule dis- 
sociating into fractionally charged species. As this 
cannot happen, the KS potential develops a step 
in the bonding region, which re-aligns the atomic 
homo's, so preventing any bias. Another way to 
put this, is that the chemical potential must be the 
same throughout the long-range molecule, and equal 
to the molecular HOMO orbital energy. Since the 



chemical potential of the true system is the small- 
est ionization potential in the system at infinite sep- 
aration, the KS potential near the atom with the 
larger atomic ionization potential must be uniformly 
raised by A/ to bring it to the ionization potential of 
the other atom, while asymptotically stepping back 
down to zero. 

We shall now introduce our two-electron model to 
study these features further and how they develop 
as a function of bond-length. 

III. A ONE-DIMENSIONAL 
TWO-ELECTRON MODEL OF LIH 

A simple one-dimensional, two-electron model of 
lithium hydride can be used to illustrate several 
important features of hetero-atomic dissociation. 
Much of the essential physics of the dissociation pro- 
cess may be captured by focusing on the chemically 
important valence electrons, while representing the 
effect of the core electrons by an average effective 
potential, such as a pseudopotential or frozen-core 
approximation. In the case of LiH, the two core elec- 
trons are localized in the Li Is shell, while the two 
valence electrons are delocalized across the molecule. 
Our goal is to analyze the effect of bond breaking 
and formation on the various xc components (Eq.4), 
which in a real molecule will be partially obscured by 
shell structure and other many-electron effects from 
the electrons in the Li Is core. Our two-electron 
model enables us to circumvent this complication, 
by focusing solely on the electrons involved in the 
bond, and their effect on the Kohn-Sham character- 
istics. As further simplication, it is reasonable to 
use a one-dimensional model, where the coordinate 
is taken to be along the bond axis, for cylindrically 
symmetric systems such as a diatomic molecule. 

As is often done in one-dimensional models, the 
Coulomb potential ±l/|r — r'| is replaced by a 
soft-Coulomb potential, ±l/-\/a -f (x — x')^. For a 
model of LiH at interatomic separation R, we write 
the electron-nuclear potential as: 



^a+{x-R/2f ^h+{x + R/2Y 

(16) 

The "softening parameters" a and 6 are directly re- 
lated to the ionization potentials of the individual 
atoms (see shortly). Similarly, the electron-electron 
repulsion is represented by a soft-Coulomb form: 

= , / == (17) 

-^C -I- [Xi - X2y 

We place Li at — i?/2 and H at R/2, and choose 
the parameters a = 0.7, b = 2.25 and c = 0.6, for 
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reasons explained in the following. With a = 0.7, 

the ionization potential of hydrogen in our model 
comes out to be 0.776H as compared with 0.5H for 
the real atom. Taking b = 2.25 yields that of Li in 
our model as 0.476H as compared with 0.198H for 
the real lithium atom. The correct difference in ion- 
ization potentials of the atoms AI = Ir—Iu = 0.3H 
is however exactly reproduced by our parameters; 
A/ is a key quantity in our analysis of the KS po- 
tential at large interatomic separations. Due to the 
long-range nature of the soft-Coulomb interaction, 
we choose the atomic ionization potentials be larger 
than in the true 3D case to prevent the atomic den- 
sities of the individual atoms from being too diffuse. 
Other factors considered were the equilibrium bond 
length (model 1.6a.u., true B.Oa.u), dissociation en- 
ergy (model 0.068 H, true 0.092H) and molecular 
first ionization potential (model 0.51H, true 0.29H), 
where the nuclear-nuclear interaction is modelled by 
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V(a + & - c) + R? 



(18) 



In Fig.(l) the dissociation curve for our model is 
plotted for comparison with that of 3D LiH^^. 
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FIG. 1: Binding energy for: 1) model ID LiH 2) true 3D 
LiH 



IV. NUMERICAL SOLUTION AND 
RESULTS 

Wc use a standard Rungc-Kutta differential equa- 
tion solver, as implemented in the octopus code'*^"'*'* 
to numerically solve for the ground-state wavefunc- 
tion ^{x,x') of the Hamiltonian: 



H 



I d2 



1 d2 



nian is mathematically equivalent to that of one par- 
ticle moving in the two dimensional potential^^"'*^: 



2da;i2 2dx2^ 



where Ucxt(a^) and Vcc{xi — X2) are defined in 
Eqs. (16) and (17). The above two particle Hamilto- 



Vext(a;) + Vey,t{y) + Vee{x - v) 



(20) 



We solve the equivalent one-particle Schrodinger 
equation on a rectangular two dimensional 25 by 25 
a.u. real space grid. The grid points are separated 
by a distance of .04 a.u. 

The density is obtained from the wavefunction 
through p{x) = 2/ dx'Yi/{x,x')\^ , and then substi- 
tuted into Eq. (12) to yield the exact KS poten- 
tial. The xc potential can be isolated from subtract- 
ing the external potential Eq. (16) and the Hartree 
potential, v-a{x) = J p{x)vcc{x — x')dx' (using Vce 
from Eq. (17)). Because Vx{x) = —Vh{x)/2, wc may 
also extract the correlation potential alone Vc{x). 
From the conditional probability amplitude 6, we 
construct fx°'°(.i;) according to Eq. (7). 

The exact KS potential is plotted at several dif- 
ferent internuclear distances in Fig. (2) alongside the 
external potential and the density. As the molecule 
dissociates, step-like and peak-like features clearly 
develop in the KS potential. There is a build-up in 
the KS potential around the more electronegative 
atom that, at each R, eventually returns to zero on 
the right-hand side of the atom (one sees the begin- 
ning of the return to zero at the smaller separations 
shown, but at separation R = 10.0 this occurs be- 
yond the region plotted). 
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(19) 



FIG. 2: Vs (solid curves), Vext, and the density (dotted 
curves) plotted at the internuclear separations indicated. 

These features occur in the response and kinetic 
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components of the correlation potential, as is evident 
in Fig. 3. Here, wc plot the xc potential (solid), 
which is the sum of the xc-hole potential (dotted) 
and the response components Wc,kin + ^'resp (dashed). 
At equilibrium bond-length {R = 1.6a.u.) the xc po- 
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FIG. 3: The total xc potential Wxc (solid curves), 
«c,kin + frcsp (dashed curves), iixc" (dotted curves) at var- 
ious internuclear separations 



tential is dominated by the potential of the xc hole. 
As the molecule dissociates the Wc,kin + i^rcsp com- 
ponents become large giving rise to clear peak and 
step structures. At large separation, the local vari- 
ation of the total xc potential around each atom is 
almost entirely due to the xc hole: at i? = 6.0a.u. 
and 10.0 a.u., v^c and v^'^'^ exactly coincide near 
Li, while near H they have the same shape, but the 
well in is translated upward by exactly 3.0 a.u. 
relative to v^"^^, which is the magnitude of the step 
A/ (see earlier Sec. II). At R = 6.0a.u., the step has 
reached its asymptotic value of Ih — Ili = 3.0a.u. 
As the molecule is pulled apart further, the step does 
not increase in size, but becomes flatter and larger 
in spatial extent. In Section IV A, wc show that 
the bond-length at which the step begins to develop 
is related to the position of the avoided crossing be- 
tween the ground-state and the state that eventually 
becomes the lowest charge-transfer state. 

A sharp peak near the rise of the step is evident 
in the xc potential (Fig. 3); this occurs in the ki- 
netic component to the correlation potential, Uc,kin, 
as discussed in Sec. II. Wc return to an analysis of 
its magnitude and location in the widely separated 



limit, and an explanation of its role in achieving the 
exact density of the interacting system, in Sec. IV B. 

In Figure 4, we plot the potentials for the sepa- 
ration R = 10.0. This is the largest separation for 
which we could converge our numerical method. In 
the limit of very large separation, we expect that the 
KS potential reduces simply to the external potential 
in the region of the nuclei, because it would be a one- 
electron system around each nucleus. There may be 
a possible shift up or down relative to the external 
potential, since constants in potentials have no phys- 
ical relevance. That is, we expect the Hartrce-plus- 
xc potential becomes flat in the atomic regions. We 
notice in our model at R=10.0, that this is approxi- 
mately true: there is however a gentle slope in tinxo 
upward around the left nucleus, and downward on 
the right, and this is largely due to a Hartree effect. 
Compared to the atomic densities in true Coulomb- 
interacting systems, the soft-Coulomb densities in 
one-dimension fall off much slower away from their 
nuclei, resulting in a longer-ranged Hartree and xc 
potential than in the true 3D counterpart. It is clear 
from the graph that the Hartree potential is still sig- 
nificant in the interatomic region. In addition to 
long-ranged correlation effects from the density on 
the "other" atom (i.e. the peak and step), the xc po- 
tential must cancel the local Hartree potential: the 
exchange potential takes care of half of this cancella- 
tion (Eq. 14), but the correlation potential must also 
contribute a well of half the size of the Hartree, as is 
evident in the graph. Despite the long-rangedness of 
the Hartree potential, R = 10.0 can still be viewed 
as "asymptotic" from the point of view of the peak 
and step structures in the correlation potential: the 
graph shows clearly that the potential on the hydro- 
gen nucleus on the right is raised by A/, and the 
peak has a height of about 0.76 (see last section). 
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FIG. 4: Components of the potentials for R = 10.0 
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A. Onset of the step: Relation to Potential 
Energy Surface Crossings 

We now show that the bond-length at which the 
step begins to significantly develop is correlated with 
the position of the avoided crossing in the potential 
energy surfaces associated with the ground state and 
the lowest excited charge transfer state. In a dia- 
batic picture, ionic and covalent curves cross at an 
internuclear distance, Rq, which is approximately 
equal to ^/{Id — ^a), where Id is the ionization 
energy of the donor and Aa the electron affinity of 
the acceptor, in the lowest charge-transfer state of 
the long-range molecule^^. When one considers the 
adiabatic potential energy surfaces, the crossing be- 
comes an avoided one, whose splitting exponentially 
decreases as a function of ii!c^^- 

What has not been previously pointed out, how- 
ever, is that the step structure in the KS potential 
begins to develop in the vicinity of the avoided cross- 
ing. Why this must be so lies in the fact that the 
step is an asymptotic feature, that arises once the 
two atoms are independent systems, and its shift of 
the eigenvalues of the more tightly bound atom en- 
sures that the ground state solution of the KS poten- 
tial has exactly half the density (i.e. one electron) 
on either side of the midpoint (see Sec. II). The 
development of the step must therefore track the in- 
dependence of the two atomic systems (measured, 
for example, by their indifference to a perturbation 
on the other atom). The avoided crossing marks 
the point at which the molecule transitions (mov- 
ing from short bond distances to longer ones) from 
a single system to two independent systems. The 
width of this transition tracks the magnitude of the 
ground-excited energy gap at the avoided crossing, 
i.e. it should be wider when the avoided crossing 
is at small bond distances and sharper when the 
avoided crossing occurs at large distances. 

Our model demonstrates this explicitly. Figure 5 
presents the ground- and first excited-state poten- 
tial energy surfaces for three different values of the 
electron-electron soft-Coulomb parameter, c. As c 
increases, the avoided crossing moves out and be- 
comes sharper; the lowest energy gap therefore de- 
creases, indicating that the transition from ionic to 
covalent character occurs more abruptly. 

In Figure 6, we plot the Hartree plus xc potential, 
w Hxc {x) = Wh {x) -f Wxc {x) for a range of internuclear 
separations i?, for c = 0.6. As this is the net po- 
tential that gets added to the external potential, we 
expect that in the limit of wide separation, it be- 
comes flat around each nucleus, since it should de- 
scribe essentially two one-electron systems. We see 
this in the graph, where a definite step is visible from 
R = 5.0 and higher. We see that it is indeed in the 
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FIG. 5: Ground-state and first-excited state (charge- 
transfer) potential energy surfaces for our model with 
c = 0.6 (top left),c = 1.0 (top right),c = 2.8 (bottom 
left). The energy differences between the surfaces are 
shown in the bottom right figure; their minimum lies at 
the avoided crossing. 



approach to the avoided crossing, at about R=4.0, 
that a shoulder first becomes clearly visible around 
the atom with the higher IP; this develops fully into 
a step of size A/, as the molecule dissociates. 




FIG. 6; The Hartree-exchange-correlation potential, 
«Hxc(2;) for our LiH model (c = 0.6); the values of inter- 
atomic separation R are indicated. 

In Figure 7, we plot Whxc for c = 2.8. The step 
begins to develop at larger i?, corresponding to the 
larger Rc where the avoided crossing occurs. Also, 
as the avoided crossing becomes sharper, the onset 
of the step happens more rapidly. 

The long-rangedness of the density in soft- 
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FIG. 7: As in Fig. 6 but with c = 2.8 



Coulomb systems means that the Hartree and 
exchange terms decay slower than in the usual 
Coulomb case. To clarify the step and peak struc- 
tures, we plot just the kinetic plus response term 
in Figure 8 for our LiH model of c = 0.6, and in 
Figure 9 for c = 2.8. The relation between the R 
at which the step develops and the avoided crossing 
discussed above is seen more clearly in these figures. 
Finally, in Fig. 10, we plot the value of Wc,kin + t^rcsp 
at the location of the atom with the larger IP (the 
H atom in our model) , as a function of internuclear 
separation i?, for various different c-values. This 
graph shows quite clearly that the development of 
the step tracks the location and sharpness of the 
avoided crossing: the larger the separation at which 
the avoided crossing occurs (i.e. larger c- value), the 
consequently larger R the step is onset, and that 
the step develops more sharply, corresponding to the 
sharper avoided crossing at larger distances. 
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FIG. 8: The kinetic and response components of the 
correlation potential Uckin+'^c.rcsp for our model with c = 
0.6; the values of interatomic separation R are indicated. 





w 

1,1 
111 








1) 
'V 












10 




p 




\7 



-15 -10 -5 5 10 15 

X 

FIG. 9: As in Fig. 8 but with c = 2.8 
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FIG. 10: The value of Wc,kin + Urosp at the location of the 
H atom in our model, as a function of the internuclear 
separation R, and with c-values as indicated. 



B. The asymptotic separation limit, and the 
significance of the peak 

Analytic expressions for the xc potential and its 
components in our two-electron model can be found 
in the separated-atom limit, by adopting the Heitler- 
London form for the wavefunction: 

^hl[x,x) = . = (21) 

where Sh.li is the overlap integral: 

Sh,Li = J (l)H{x)4>Li{x)dx (22) 

We will focus on the interatomic region, far from 
either nuclei, in this limit. To lowest order in the 
separation R, the orbitals 4>Li{x) and 4)h{x) in this 
region may be written: 

^Ux) = V^e-""("+^) (23) 
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where a = v2/, with I being the first ionization po- 
tential of the atom. Similar expressions hold for the 
three-dimensional case, with Coulomb interaction; 
the only differences being that instead the orbitals 
have asymptotic dependence according to 



(24) 



(where the x-axis is taken to be the bond axis). 

It is a simple exercise to construct the first-order 
density-matrix and the density using these orbitals. 
Substituting into Eqs. (12) and (9) yields the large- 
separation limit of the KB potential and VkinCi")- 

In the limit of large interatomic separation, the 
Hartree potential vanishes as the inverse distance 
from the nuclei in the inter-atomic region. Also, in 
this limit, the second order density matrix factorizes 
into a product of densities and it follows from Eq. (5) 
that f xc (t) also falls off as the inverse distance from 
the nuclei in the interatomic region. The KS poten- 
tial is then dominated by contributions from Vkin{r) 
and Wrosp(?")- Explicitly, in one dimension this is 
given by: 



' |2 



2 |(/'H|2 + |0L»P + 2Vi(/'ff0z 



2 \(l)H? + \<l>Li? + '^\r^MLi 

4 {\<j>H? + \<l>Li?+2^e<^H<t>Li? 



In the above expression 0' and (f)" dc;note the first 
and second spatial derivatives of the orbital, and e 
is the square of the overlap integral at interatomic 
separation R: 



e = 



OlHOlLi 



[an - OLLiY 



-RotLi 



— e 



(26) 



In Fig. (11), the asymptotic expression for Uhxc(= 
^^s — ^^ext) using the orbitals of Eqs. 23 is plotted for 
comparison with the Wc,kin('') + Vresp{r) component 
of the numerical solution using the soft-Coulomb po- 
tentials. (As noted earlier, the soft-Coulomb orbitals 
are longer-ranged than their 3D Coulomb counter- 
parts, so Whxc achieves its asymptotic form only at 
larger distances.) We see that the step reaches its 
asymptotic limit more quickly than the peak. For 
instance, at i? = 10.0 a.u. the peak for the numer- 
ical solution is somewhat smaller than that of the 
analytic expression, although the step has already 
reached its asymptotic value of 3. OH. 

We next derive asymptotic expressions for the lo- 
cation and magnitude of the peak and step struc- 
tures as functions of the internuclear separation. 



Defining the location of the peak from the condition 
0, we obtain 



Cpeafc 



^peak — 




(27) 



where In and Ih are respectively the ionization po- 
tentials of Li and H. For the 3D case the second term 
on the right is modified to be: 



Jh 



(28) 



Defining the location of the step by its inflection 
point, i.e. 

tains the same result, i.e. 



j2 

from the condition ^ 



"^^resp = 0; One ob- 



^step — ^peak 



(29) 



Therefore, in the asymptotic limit, our two-electron 
model shows that the location of the peak and step 
coincide. The second term in Eq. (27) is negative, 
but in general small compared to the first term for 
large inter-atomic separation R. Therefore, the peak 
and step structures are located closer to the hydro- 
gen atom; more generally, closer to the more elec- 
tronegative atom of a diatomic molecule. On the 
other hand, the minimum of the density: 



.(n) 




(30) 



lies closer to Li than the peak/step location, but still 
on the side of the bond mid-point closer to H: The 
first term of Eq.30 is identical to Eq. (27), while the 

second term contains the logarithm of the ratio -j^ 

instead of the ratio 4^ , which is smaller than one. 

Our simple two-electron model thus explains the 
earlier observations in real molecules^^: In the gen- 
eral many-electron hetero-atomic case, given that 
the peak and step structures arise from the delocal- 
ized HOMO, our analysis can predict their positions. 
The location of the step was seen to coincide with 
the peak in the true LiH molecule, with both lying 
closer to the H atom, at least for the largest inter- 
atomic distances that those calculations were able 
to perform. For the homo-atomic case, our results 
(Eqs. 27 and 30) predict that Xpea.k = a^dens.min = 
and so the minimum of the density and peak loca- 
tion coincide at the bond midpoint; also borne out 
by the examples in the literature. 

We next turn to the magnitudes of the structures. 
Using the density matrix constructed from the or- 
bitals in Eqs. 23, one can show that the magnitude 
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of the peak structure in Vkin{r), in the Umit that the 
overlap integral vanishes, is given by the expression: 

<kT„ = J(v^+/^)' (31) 

For our two-electron model of LiH, this gives a value 
of 0.616 a.u. Adding the value of the step in Vresp 
at its inflection point (A//2 = 0.15au), gives .7672 
a.u., which is indeed what the peak of our numerical 
solution asymptotes to. For the homo-atomic case, 
the above expression gives a value of V^eak = 0.5au, 
agreeing with the results of Rcfs.'^^ and-'^^ for the 
true homo-atomic two electron system H2. How- 
ever, in Ref.^^, the magnitude of the peak for true 
LiH, was significantly smaller than this prediction. 
This discrepancy is due to the effect of the localized 
core electrons in the Li Is shell, which lead to a dra- 
matic decrease in the magnitude of the gradient of 
the conditional probability amplitude eq. (6) in the 
inter-atomic region, and hence by eq.( 9), a decrease 
in the magnitude of the peak. 
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FIG. 11: Asymptotic expression for Uhxc (solid curve) in 
the inter-atomic region compared to Vhxc in our model 
(dashed curve), and Ukin+Wresp (dotted line) in our model 
at separations indicated. 

As discussed in Sec. II, the peak emerges out of 
analyzing the change in the conditional probability. 
We now give a different argument for why the peak 
must be there, even though it has negligible effect on 
the ground-state energetics. The peak occurs when 
one takes the "non-bonding" orbital as the KS or- 



bital: 

</> = ^/{pH + pLi)/2. (32) 

(Here Ph is the atomic density of the H atom and pLi 
that of the Li atom, i.e. the squares of the orbitals 
in Eq. 23). This is the exact doubly-occupied KS 
orbital, since twice its square yields the exact density 
in the limit of infinite separation, p = 2\(j)\'^ = pjj + 

PLi- 

If one instead takes the " bonding orbital" : 

<^bond = (VpW2 + VpLi/2), (33) 

and finds the KS potential corresponding to this, 
there is no peak structure (but there is still the step). 
That is, if one asks what is the KS orbital for the 
KS potential with the peak structure sliced out, the 
KS orbital would instead be 4>hond- Now the density 
corresponding to ^!>bond is 

Pbond = 2|(/)bond|^ = pH + PLi + "^^/PhPLi (34) 

i.e. is equal to the sum of the atomic densities 
plus a term 2^pHPLi- This term is indeed very 
small, but taken as a fraction of the total density, 
^JPHPLi/ (ph + PLi)i displays a peak at the exact 
same location as the peak in the exact KS poten- 
tial, Eq. 27 (Figure 12). The shape of the peak is 
different but its maximum coincides in the limit of 
infinite separation. This suggests an interpretation 
of the peak in tic, kin (in the exact KS potential), as 
a barrier that pushes back to the atomic regions the 
extraneous density 2^pHPLi that would be in the 
bonding region if the peak was absent. Since the KS 
system by definition must get the density correct the 
peak must be there. 

The interpretation here is closely related to the 
analysis of Ref.^^ of ho mo- atomic molecules, where 
it was shown that the kinetic energy density for the 
exact KS orbital develops a well in the bond mid- 
point region, that must be compensated by a peak in 
the KS potential in order to keep the constant value 
of the KS orbital energy. An LCAO approximation 
to the orbital (analogous to (?ibond above) does not 
display the well. 

V. DISCUSSION AND IMPLICATIONS 

Using a simple one-dimensional model of a two- 
electron heteroatomic molecule we studied features 
of the exact KS potential that arise for real 3D het- 
eroatomic molecules. In particular we examined the 
characteristic step and peak structure in the inter- 
nuclear region, that develop as the molecule dis- 
sociates. These unusual features are a peculiarity 
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FIG. 12: The peeik in Vc,kin{x) and that in the fractional 
density-error of the KS orbital solution to the KS poten- 
tial with the peak taken out, y/pHPLi/{pH + Pu), (see 
text) have the same location. 

of the non-interacting KS description: on the one 

hand, as a molecule dissociates, the interaction be- 
tween the electrons on one atom and those on the 
other vanishes, so why, fundamentally, do such stark 
structures appear in the KS potential? The answer 
ultimately lies in the single-Slater-determinant de- 
scription in the KS system: although this is indeed 
how the exact KS system describes the state, it is 
far from the true wavefunction which needs, even 
qualitatively, two Slater determinants. In the two- 
electron model, the KS system consists of a doubly- 
occupied spatial orbital, blatently far from the true 
two-orbital interacting system. Mathematically, the 
structures can be understood by considering the re- 
sponse and kinetic components of the correlation po- 
tential, as explained in earlier works and in Section II 
of the present paper. Physically, a KS potential that 
lacks the step leads to dissociation into fractional 
charges; a KS potential that lacks the peak leads 
to a KS orbital that yields an incorrect (albeit ex- 
ponentially small) density in the intcrnuclcar region. 
The former point is well-recognized in the literature, 
while the latter point elaborates on an earlier inter- 
pretation^i (Sec. IV B). 

Due to the simplicity of our two-electron model, 
we are able to investigate in much more detail than 
in the earlier literature, the development of these 
structures and their asymptotic properties. Several 
of these features carry over to the true many-electron 
3D case, since they arise from the HOMO orbital. 
We showed that the step begins to develop at the 
internuclear separation where the avoided crossing 
in the ground and lowest charge-transfer state is ap- 
proached, and explained why. We gave an exact for- 
mula for the location of the step and peak, in the 
limit of large separation, finding that the two struc- 
tures are located at the same place, and closer to 



the atom with the larger IP, consistent with the few 
calculations done on real molecules in the literature. 

Being in a region of very low electron-density, 
these features, in themselves, have little energetic 
consequences for the ground states of these sys- 
tems. However they have dramatic consequences 
for time-dependent processes, excitations, and re- 
sponse. For example, it has been shown that the 
related peaks that appear in the interatomic regions 
of a hydrogen chain significantly (and correctly) re- 
duce the polarizability of the chain and that local 
and semi-local approximations which lack the peak, 
consequently significantly underestimate the (hyper- 
) polarizability^^. As TDDFT begins to be utilized in 
molecular transport calculations, we anticipate the 
peaks will act as barriers decreasing the current. 

The step in the KS potential ultimately imposes a 
rather complicated structure on the exact xc ker- 
nel of TDDFT37.38^ Because of the reah gnment 
of the atomic HOMO's, the molecular HOMO and 
LUMO are symmetric and antisymmetric combina- 
tions of the atomic HOMO's, separated in energy 
merely by the tunnelling factor, that vanishes as 
exp(— const. i?) as the molecule dissociates. There- 
fore three KS determinants become near-degenerate: 
the doubly-occupied HOMO, a single-excitation to 
the LUMO, and a double-excitation to the LUMO. 
That is, the step introduces static correlation in 
the KS system that is not present in the true in- 
teracting system. It is the job of the TDDFT xc 
kernel to "undo" this static correlation, in order 
to yield good excitation energies in the true sys- 
tem. This has a dramatic effect on the structtirc 
of the xc kernel for charge-transfer excitations in 
molecules composed of open-shell fragments^'''^^; in 
particular, the double-excitation induces a strong- 
frequency-dependence on the kernel. 

Almost all the approximations in use today do not 

capture the step and peak structure in the potential. 
Carefully constructed orbital functionals for the cor- 
relation potential may display these structures, as 
has been explicitly shown in Ref.''®. Interestingly, 
static correlation in the KS system is nonetheless not 
escaped in the usual (semi-)local approximations. 
Delocalizcd orbitals underlie the fractional charge 
problem, and the HOMO and LUMO become near- 
degenerate as the molecule dissociates. Figure 13 
demonstrates this for the LiH molecule within LSD; 
a similar merging of the HOMO and LUMO is also 
seen in GGA. 
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FIG. 13: LSD KS eigenvalues for LiH become near- 
degenerate as a function of internuclear separation R 
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